shp <- readOGR(paste0(droot,
"data/01 - Getting Started/vg250_2019-01-01.gk3.shape.ebenen/vg250_ebenen/VG250_KRS.shp"))## OGR data source with driver: ESRI Shapefile
## Source: "C:\Users\besel\OneDrive\Academia\Research\TLS\02 - First Draft\data\01 - Getting Started\vg250_2019-01-01.gk3.shape.ebenen\vg250_ebenen\VG250_KRS.shp", layer: "VG250_KRS"
## with 431 features
## It has 23 fieldskreise.shp <- sp::spTransform(shp, CRS=CRS("+init=epsg:4839"))
kreise.df <- fortify(kreise.shp, region = "AGS")empl.counts <- read_csv(paste0(droot, "output/01-First Data Exploration/empl_count_DDRWZ_KR2019.csv"))
empl.counts2 <- read_csv(paste0(droot, "output/01-First Data Exploration/empl_count2_DDRWZ_KR2019.csv"))
crswlk <- read_dta(paste0(droot, "data/01 - Getting Started/Crosswalk Counties/kr_1989_ost_2019_weight.dta"))
mat <- empl.counts %>%
column_to_rownames(var = "kr_201901") %>%
replace(is.na(.),0)
eci1 <- KCI(mat, RCA = TRUE)
eci1 <- data.frame(empl.counts$kr_201901, eci1)
colnames(eci1)[1] <- "id"
eci1 <- eci1 %>%
left_join(y= crswlk[, c("kr_201901", "name_2019")], by=c("id" = "kr_201901")) %>%
mutate(id=replace(id, id==11200, 11000)) # make the Berlin change
tt <- kreise.df %>%
mutate(id = as.integer(as.character(id))) %>%
left_join(y=eci1, by="id")
# Merge with total employment per 2019 county
tt <-tt %>%
left_join(y=empl.counts2[, c("kr_201901", "totsum")], by=c("id" = "kr_201901"))g <- tt %>%
#filter(eci1!=is.na(eci1)) %>%
ggplot(aes(x=long,y=lat, group=group, fill=eci1, name=name_2019, info=totsum))+
geom_polygon()+
ggthemes::theme_map()
#coord_map()And then do plotly (plotly file is to big to show in html. has to be exported)
# Get correlation between totsum and eci
sk <- tt %>%
distinct(eci1, .keep_all = TRUE) %>%
drop_na()
cr <- round(cor(sk$totsum, sk$eci1), 2)
# Make interactive plot
ggplotly(g, tooltip = c("name_2019", "eci1", "totsum")) %>%
layout(title = list(text = paste0('1989 ECI on 2019 county map',
'<br>',
'<sup>',
'Correlation between ECI and total employment is ', cr,
'</sup>')))